IMPLICIT METHODS FOR EQUATION-FREE ANALYSIS: 
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Abstract. First, we give a rigorous convergence result for equation-free analysis in the setting 
of slow-fast systems using implicit lifting. Second, we apply this result to study the idealized traffic 
modeling problem of phantom jams generated by cars with uniform behavior on a circular road. It 
is shown, that the implicitly defined coarse-level time stepper converges to the true dynamics on the 
slow manifold with an accuracy that is beyond all orders of the small parameter measuring time scale 
separation. These results are applied to investigate the behavior of the microscopic traffic model on 
a macroscopic level. The traffic jams are waves that travel slowly and in opposite direction compared 
to the car velocity. The standard deviation is chosen as a macroscopic measure of traveling wave 
solutions and is continued on the macroscopic level in the equation-free setup. The collapse of the 
traffic jam to the free flow corresponds in the relevant parameter region at the macroscopic level to a 
saddle-node bifurcation of the traveling wave. We continue this bifurcation point in two parameters 
using equation-free analysis. 
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1. Introduction. When one studies systems with many degrees of freedom, for 
example, systems with a large number of particles or interacting agents, one is often 
interested not so much in the trajectories at the microscopic level (that is, of individ- 
ual particles) , but in the behavior on the macroscopic scale (of the overall distribution 
of particles). The classical example are the motions of molecules of a gas, resulting in 
the laws of thermodynamics. In this classical case the macroscopic description is de- 
rived in statistical mechanics from knowledge about the microscopic behavior through 
time scale separation. Other important examples are emerging patterns in physical, 
chemical and biological systems, e.g., Rayleigh-Benard convection rolls [35], Belousov- 
Zhabotinsky reaction [3, 44J , stripes on zebra skin or pattern on butterfly wings |41j . 
A common approach in physics literature to deriving macroscopic descriptions are 
the so-called adiabatic elimination or the slaving principle [HI [16] . These concepts 
are related to the theorems about center-manifold and slow-manifold reduction of the 
mathematical literature [7[ HOI HO] • 

For systems where no explicit macroscopic description can be derived from micro- 
scopic models, Kevrekidis et al proposed that, if the number of particles is moderate, 
then it is sometimes possible to skip the derivation of a macroscopic description by 
performing the analysis of the dynamics in the macroscopic scale directly. This ap- 
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Fig. 1.1. Sketch of the macroscopic time stepper <J>(t; ■) . TTie macroscopic state x (to) is mapped 
to a microscopic state u(to) by using the lifting operator C. The available microscopic time stepper 
is used to evolve the system to the microscopic state «(to + t), which is mapped to a macroscopic 
state x(tg + t) using the restriction operator 1Z. This procedure constitutes the coarse-level time 
stepper &(t; ■). 



proach relies on evaluating short bursts of appropriately initialized simulations of the 
microscopic model (see, for example, (HJ [331 HI] f° r recent reviews). It is called 
equation-free, because it assumes that the macroscopic model exists, but is not avail- 
able as an explicit formula. Equation-free methods are particularly appealing if either 
explicit macroscopic descriptions are unavailable, or one wants to study the underly- 
ing system near the boundary of validity of its macroscopic description (for example, 
as one decreases the number of particles, finite size effects may start to appear as 
small corrections to the macroscopic model) . Equation-free analysis has been applied 
for a large class of multi-scale models that fit roughly the description of singularly 
perturbed systems [9] in a broad sense (see motivation in [24] ) . such as stochastic 
systems (28] |37] , agent-based models [5j El [14] , molecular dynamics [4] or neural dy- 
namics [26, 33 , to perform high-level tasks such as bifurcation analysis, optimization 
or control design [SI [3S] . 

The basic building block of equation-free analysis is an approximate coarse-level 
time stepper &(t; •) for short times t (compared to the slow time scale) in the phase 
space of macroscopic variables (say, K rf ). This coarse-level time stepper is typically 
defi ned with the help of three operations: lift C, evolve and restrict 1Z (cf. Figure 



1.1). To compute the map $>(t;x) on a given macroscopic state x £ M. d , one has to 
apply a lifting operator C to map x to a microscopic state u € M. D (typically, D 3> d), 
then one runs the microscopic simulation for the time t, and finally one maps the 
end state of the microscopic simulation back into M. d using a restriction operator 1Z. 
A proof of any claim that this would be a good approximation of the true dynamics 
of the macroscopic variable x in a given example will have to invoke the following 
sequence of arguments. Initially assume that the microscopic system is a slow-fast 
system with a transversally stable slow manifold, for which the macroscopic quantity 
x is a coordinate. The first question is then: does the approximate coarse-level time 
stepper $ converge to the true dynamics on the slow manifold in the limit e — » 0, 
where s is the parameter measuring the time scale separation? In addition to the case 
discussed here, equation-free analysis is also applied to high-dimensional, stochastic 
(or chaotic) systems showing macroscopic behavior because the microscopic degrees 
of freedom 'average out rapidly' [H GUI El| ■ I n these cases another question must be 
addressed: in which sense is the averaging process approximating a classical slow-fast 
system? 

1.1. An implicit coarse-level time stepper. A prerequisite before equation- 
free analysis can be performed is the problem of finding the restriction and lifting 
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operators TZ and C. Figure [O] suggests the relation <E>(i; •) = TZ o evolve o £, which 
is not justified in general. Let us assume that the microscopic system is slow-fast 
and the macroscopic system corresponds to the slow flow on the slow manifold in the 
coordinate x. Then it is clear that an arbitrary choice of C and TZ does not lead to a 
coarse time-stepper $ which approximates the slow flow in any way, even in the limit 
of infinite time-scale separation (e — > 0). The source of the error is an initialization of 
the microscopic system away from the slow manifold. One relies on the separation of 
time scales in a so-called healing step to reduce this error. However, in most reviews 
this healing is applied inconsistently [HI 1231 HI] (that is, healing would not lead to 
$ converging to the true slow flow in the limit of infinite time-scale separation, even 
in the ideal case of a slow-fast system). A consistent way to perform healing are 
so-called constrained-runs corrections after lifting, developed in [T5J 22J 23] . These 
papers developed schemes of increasing complexity to compensate for this error source. 

An alternative, explained in Section [2] is to use an implicitly defined coarse- level 
time stepper <E>, where the slow flow is not measured at pre-determined points in 
space, but rather at healed points. In the special case of computation of equilibria 
the use of the implicit time stepper reduces to the formula introduced as the "third 
method" by Vandekerckhove et al [35] • In Section [3j we give a detailed proof of the 
convergence of the implicitly defined coarse-level time stepper $ to the flow on the 
slow manifold, answering the question of convergence for the implicit time stepper. 
The approximation accuracy of $ (under some transversality conditions) is beyond 
all orders in the parameter e measuring the time scale separation. In Section [4] we 
discuss the assumptions and consequences of the convergence theorem and compare 
it to other results in the literature. 

1.2. Macroscopic behavior of a microscopic traffic models. In Section [5] 
and Section [6] we apply the convergence results to study a traffic modeling problem 
that fits perfectly into the framework of equation-free analysis: a large number of cars 
(the microscopic particles) on a circular road that interact with each other, resulting 
in so-called phantom jams moving slowly along the road against the direction of traffic, 
i.e., forming a traveling wave at the macroscopic level. 

The mathematical modelling and analysis of traffic flow dynamics has a consider- 
able history (see, e.g., [T7J[251|3T] for reviews). Macroscopic traffic models use partial 
differential equations, such as Burger's equation [2"5] , for modelling the flow. They 
model the density of cars as a continuous quantity to formulate directly macroscopic 
equations for density and flux along the road. In contrast, microscopic particle mod- 
els for the cars (deterministic [I] or stochastic [TH1I3S]) can be used to describe the 
behaviour of individual cars or drivers. An advantage of microscopic models is that 
parameters can be assigned directly to the individual drivers' behavior (for example, 
aggressiveness, inertia or reaction delay) such that these parameters' influence and the 
trajectories of individual cars can be investigated. Another use of microscopic models 
is to test the effect of new devices for individual cars, for example, cruise control, on 
the overall traffic prior to their implementation in real traffic. In this paper we use 
the optimal velocity model [T] as an example of an underlying microscopic model. 
The optimal velocity model results in a set of coupled ordinary differential equations, 
but despite its simplicity it can reproduce the phenomenon of phantom traffic jams. 
An advantage of choosing the optimal velocity model is that we have guidance from 
the results of direct bifurcation analysis of the full microscopic system when only a 
few cars are involved [111 130) as well as from perturbation analysis based on discrete 
modified Korteweg-de-Vries equation [10] . Direct bifurcation analysis of the micro- 
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scopic system becomes infeasible when the number of cars gets large. Furthermore, 
it is difficult to analyze macroscopic quantities for which typically no equations are 
explicitly given such as the mean and standard deviation of headways or densities of 
cars. In Section [6] we show how this difficulty can be tackled by using equation-free 
methods for the bifurcation analysis on a macroscopic level. 

In Section [7] we summarize the obtained results and give an outlook on open 
problems. 

2. Non-technical description of general equation-free analysis with im- 
plicit lifting. Equation-free analysis as described by [23J is motivated by ideas from 
the analysis of slow-fast systems: one assumes that on a long time scale the dynamics 
is determined by only a few state variables and the other state variables are slaved. 
Mathematically this means that the flow of a high-dimensional system under study 
converges rapidly onto a low-dimensional manifold on which the system is governed 
by an ordinary differential equation (ODE). In many practical applications the con- 
vergence is achieved only in the sense of statistical mechanics (the effects of many 
particles averaging out; see 0|6]). The traffic problem does not require this notion of 
weak (averaged) convergence such that we give our description and subsequent con- 
vergence proofs of equation-free analysis using terminology of slow-fast systems with 
transversally stable slow manifolds following the notation of [§]• 

2.1. The notion of a slow- fast system. Let 

u = f e (u) (2.1) 

be a smooth dynamical system defined for u € M. D , where f £ depends smoothly on 
the parameter e. We assume that e is a singular perturbation parameter. This means 



that the flow A/ e generated by (2.1 1, 



M E : K x R D -> R D , (t;u) i-> M e {t;u) 

has a whole smooth <i-dimensional submanifold Co of equilibria for e = 0: if u G Co 
then M (t;u) = u (and, thus, fo{u) = 0) for all t. The dimension d is the number of 
slow variables. In the notation of singular perturbation theory, t measures the time 
on the fast time scale. We assume that this manifold Co is transversally uniformly 
exponentially stable for e = 0, which corresponds to the stable case of Fenichel's 
geometric singular perturbation theory [9j. For this case we know that the flow 
M s (t; •) has a transversally stable invariant manifold C £ for small non-zero e, too (in 
Section [6j e ~ 1/N). This manifold C e is called the slow manifold, and the flow M e , 
restricted to C £ , is called the slow flow. 

2.2. Lifting, restriction and time stepping. The equation-free approach to 
coarse graining [23] does not require direct access to the right-hand side f e of the 



microscopic system (2.1| but merely the ability to evaluate M e (t; u) for finite positive 
times t (t <C 1/e in the fast time scale t) and arbitrary u. It also relies on two smooth 
maps that have to be chosen beforehand 

1Z : K D — > R d the restriction operator, 

C : R d -> R D the lifting operator. 

In the optimal velocity (OV) model discussed in Section [6j 1Z is chosen as a mapping 
from headway profiles to the standard deviation a and C constructs a headway profile 



by using a (cf. ( |6.3| and (6.4 1). 
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Fig. 2.1. Sketch showing a typical geometry in a slow-fast system with a slow manifold C E 
and an arbitrary lifting C and restriction Ti. The healing iWe(t s kipi') * s applied to all points in 
the domain of C. Note that dom£ and rgTZ can be different, but must have the same dimension. 
(rgTZ) 1 ^ refers to an arbitrary transversal complement of rgTZ. 



The basic idea underlying [53] is that one can analyze the dynamics of (2.1 1 on 



the slow manifold C e by studying a map in the space of restri cted variables x in the 



domain of £ (called dom£ C K d ) of the form (cf. Figure 1.1 ) 

Lift — > Evolve — > Restrict, 

or, to be precise, the map 

P E (t;-):x^Tl(M e (f,jC(x))) = [TZ o M e (t; •) o £](z) (2.2) 

for selected times t <C 1/e. The central question is: how can one compose a macro- 
scopic time stepper, that is, an approximate time-c) map $((5, •) : K d — > K d using 
coordinates in the domain of C for the flow M e restricted to C e l One important ob- 
servation is that this map <!> must be defined implicitly. Figure [2~T] shows how one can 
define a good approximate time-5 map $(<5; •). It contains an additional parameter 
iskip, called the healing time in [23 . This healing time must be applied to both, the 
argument x and the result y of Thus, $((5; •) is given implicitly by solving 

P £ (t skip ;y) =P e (i skip +S;x), that is, 

K(M £ (t skip ; £(y))) =K(M s (t skip + 6; C(x))) for y, and setting (2.3) 
$(8',x) :=y. 

Under some genericity conditions on TZ and C and M e the order of approximation for 
$ is exponentially accurate for increasing t sk i p and 6 > 0: 

\\${5;x)-$*(5;x)\\ < Cexp(-K t skip ). (2.4) 

In this estimate K > and C > are constants independent of x, S and < s ki P for x in 
a bounded open set dom£ C M. d . The flow is the exact flow M e , restricted to the 
slow manifold C £ , in a suitable coordinate representation in dom£. 

The details of the convergence statement and its conditions will be made more 
precise in Section [3j The convergence statement holds not only for the values of 
the approximate flow map $ but also for its derivatives with respect to the initial 



value. So, d^^{§;x) — d J 2 ^*{5]x) satisfies inequality (2.4), too (possibly with other 
constants C) for derivative orders j less than a given k (the subscript of d refers 
to the argument of $ with respect to which the derivative is taken). The degree of 
achievable differentiability is determined by the timescale separation: the smaller e, 
the smoother the slow manifold C e , and, thus, the higher we can choose the maximal 
derivative order k. 

Based on the implicitly-defined approximate flow map <&, one can now perform 
higher-level tasks in equation-free analysis. 

2.3. Bifurcation analysis of macroscopic equilibria. Bifurcation analysis 
for equilibria boils down to finding fixed points and their stability and bifurcations 
for $(<5; •) with some small, arbitrary 5 (that is, 6 <C 1/e in our notation). In terms of 
1Z and £, the equation <fr(5;xo) = x , defining the equilibrium x , reads (cf. Figure 
2lb 



TZ(M £ (t skip + S;£(x ))) = H(M s (t skip ;£(x ))). (2.5) 

This equation has been proposed and studied already in [35]. In applications, ( |2.5[ ) is 
solved using a Newton iteration (cf. ( |6.9| in the OV model). Since the time stepper 
is defined implicitly, one finds the stability and bifurcations of an equilibrium xo by 
studying the generalized eigenvalue problem 
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[R(M e {t skip ;£{x)))] 



(2.6) 



This eigenvalue problem will give the eigenvalues of the implicitly-known flow <I>(<5; •) 
linearized in the equilibrium xq such that bifurcations occur when A is on the unit 
circle. 

2.4. Projective integration. In projective integration one approximates the 
ODE for the flow on the slow manifold C e in the coordinate x g M. d . The ODE for the 
true flow on the slow manifold is an implicit ODE with the solution x(t), which 
will be derived in detail in Section [3] Its approximation based on $ is 



-n(M e (t skip ;£(x))) 



85' 



K(M e (t skip + 5; £(x))) 



6=0 



(2.7) 



For fixed t sk i p the left-hand side is a function of x G R d such that the time-derivative 
of this function defines (implicitly) the time derivative of x. The term inside the 
partial derivative on the right-hand side is a function of two arguments, S and x, for 
which one takes the partial derivative with respect to its first argument 5 in 5 = 0, 
making also the right-hand side a function of x only. Consequently, every integration 
scheme becomes implicit. For example, if one wants to perform an explicit Euler step 
of stepsize h starting from Xj at time tj, this is equivalent to the implicit scheme 
(defining Xj+\ as the new value at time tj+\ = tj + h): 
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(2.8) 



or, in terms of restricting and lifting, 

n{M e (t sMp ;£(x J+1 ))) - n(M e (t skip ;£(x 3 ))) 



= - [ft(M E (i 8kip + 5; C( Xj ))) - K(M £ (t skip ; £(x 3 )))} 
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Projective integration becomes attractive if one can choose h either much larger than 
t s kip and 5, or one can set h negative enabling integration backward in time on the slow 
manifold (cf. ( |6.14 ) and Figure 6.1), even though the original system is very stiff in 



Ji D (and thus, strongly expanding backward in time). Note that during computation 
of residuals and Jacobian matrices one can evaluate P e (t s kipl x) as a by-product of the 
evaluation of -P e (i s kip +8;x), assuming that the restriction 1Z is of comparatively low 
computational cost. 

2.5. Matching the restriction. Sometimes it is of interest to find a microscopic 
state u € M. D on the slow manifold C E that has a particular x £ K d as its restriction 
(H(u) = x), see [T^l H3] - This state u is defined implicitly, and can be found by 
solving the d-dimensional nonlinear equation 

TZ{M e {t skip ;C(x)))=x (2.9) 

for x, and then setting u — M £ (t sk i p ; C(x)). This solution u is close to the true 
slow manifold C e with an error of order exp(— Kt sk i p ), where the decay rate K > 
is independent of e and i s kip- This implies that, if we choose i s ki P = 0(e~ p ) with 



p G (0, 1), the distance of u to C e is small beyond all orders of e. Equation (2.9) was 
also proposed and studied in [3j3 (called InitMan in [59"]). 

3. Convergence of equation-free analysis. This section gives a detailed dis- 
cussion of the convergence results of the methods sketched in Section [2| Sections[5] 
and [6] study the OV model for traffic flow as an example of the results in Section [2j 

We formulate all assumptions on 1Z, C and M e for the singular perturbation 
parameter e at e = 0, even though it is typically difficult to vary e in complex 
model simulations. However, stating the conditions at e = guarantees that they 
are uniformly satisfied for all sufficiently small e, which is the range of parameters for 
which the statements of this section are valid (cf. [S]). Throughout this section various 
constants will appear in front of exponentially growing or decaying quantities. As the 
concrete values of these constants do not play a role we will use the same variable 
name C on all occasions without meaning them to be the same. We will state which 
quantities the constant C depends on whenever we use exponential estimates. 

3.1. Existence of transversally stable slow manifold. In order to avoid the 
discussion what happens when the flow M s reaches the boundary of the slow manifold 
C e while following the slow dynamics, we assume that the manifold Cq of equilibria of 
Mo is compact. Then the flow M e , restricted to C e , will be defined for all times t £ K. 

Our first assumption guarantees transversal stability of Co- 

Assumption 1 (Separation of time scales and transversal stability) . The flow M 
converges to the slow manifold Co with a uniform rate Kq from all initial conditions 
u in some neighborhood of Cq . That is, for every u in an open neighborhood IA of the 
slow manifold Cq there exists a point p G Cq such that 

lim Mo(t; u) = p 

{note that for e = all points on the slow manifold Cq are equilibria), and the distance 
can be bounded via 

\\M (t;u)~p\\ <Cexp(-K t)\\u-p\\, \\d° 2 M (t]u)\\ <Cexp{-K a t) 



for all t > and j > 1, where the constant C depends only on the derivative order j . 
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(a) £ = (b) £ = (c) £ > 

Fig. 3.1. Sketch of geometrical interpretation of transversality assumptions. Note that (rglZ)^ 
refers to an arbitrary transversal complement of YglZ. Panels (a) and (b) show the geometry at 
e = 0: (a) the trajectory from C(x) must converge toCo, its limit is called go(C(x)) . The overall map 
1Z o go ° C must be a local diffeomorphism. This entails that &R. must have full rank on the tangent 
space Afo(uo) in any uo 6 Co (also shown in (a)), and that rgC intersects each fiber (the set of 
points u converging to the same uq £ Co) transverally (shown in (b) for d = 1 and D = 3). Shown in 
panel (c): g s andC e are 0(e) perturbations of go andCo, and M £ (t s ki p ; C(x)) — Af £ (i s j 5 i p ; g e (C(x))) 
are exponentially close for t s kip 

> 0. 



Since the slow manifold Co is compact, one can choose a uniform constant C for 
all u in the neighborhood IA. The above assumption implies the existence of a smooth 
map (called the stable fiber projection) 

g :U -> C , defined by g (u) := p, (3.1) 



assigning to each u its limit p G Cq under the flow M (see Figure 3.1 ^a)). 

We recall now two central persistence results of classical singular perturbation 
theory [S]. First, the slow manifold Cq persists for sufficiently small e, deforming to 
a smooth nearby manifold C £ (as shown in Figure [3~T| c)). Restricted to C £ the flow 
is governed by a smooth ODE (the slow flow) with a right-hand side for which all 
derivatives up to a given order k are proportional to e (larger k requires smaller e): 

\\fe(u)\\ < Ce, \\d j f E (u)[v 1 ,...,v j ]\\<Ce\\v 1 \\-...-\\v j \\ (3.2) 

for all j = {1, . . . , k}, u £ C e and Vi,...,Vj € A/" e (u) (J\f e (u) is the tangent space, 
for e = it is the null space of the linearization, of the slow manifold C £ in u). The 
constants C are independent of u. Thus, the flow M e (t; •) is a global diffeomorphism 
on the slow manifold C £ which has growth bounds of order e forward and backward 
in time: 

||^'M e (t;-)kll < Cexp(e\t\), H^Mf^f; -)k II < Cexp(e\t\), (3.3) 

for some constant C independent of t and e and all derivative orders j up to a fixed 
order k. Note that M j r 1 (i; •) = M £ (— t; ■) for all times t as long as one restricts the 
flow M e to the slow manifold C e . 

Second, the stable fiber projection map go persists for small e, getting perturbed 
smoothly to a map g e , defined for each u in the neighborhood U of the slow manifold 



Co (and its perturbation C e ). The map g £ picks for every point u ElA the unique point 
g e {u) inside the slow manifold C £ such that the trajectories starting from u and g e (u) 
converge to each other forward in time with an exponential rate K of order 1 (that 
is, K is uniformly positive for all sufficiently small e and all u € IA): 

\\4M e (t;u) - 4M e (t;g e ( U ))\\ < C exp(-Kt)\\u - g e (u)\\, (3.4) 

for all t > 0, u € U and < j < k, where the constant C is uniform for u G U. 
In general, the decay rate K has to be slightly smaller than the rate Kq asserted to 
exist in Assumption [I] for e = 0. More precisely, for every rate K < Kq there exists a 



range (0,£o) of e for which (3.4) holds. Choosing Eo smaller permits one to choose K 
closer to K . In ( |3.4| the notation d 3 2 M e refers to the jth-order partial derivative of 
the flow M e with respect to its second argument (the starting point), and the zeroth 
derivative refers to the value of flow M s (t; •) itself. The stable fiber projection map 
g £ is an order-e perturbation of go: 

\\&g E (u)-cPg (u)\\ <Ce (3.5) 

for all j = {0, . . . , k} and a constant C that is uniform for all u G U. The black curves 
transversal to C e in Figure |3.1[ c) illustrate the fibers, that is, which points of U get 
mapped onto the same point in C £ under g e . Note that the fibers are not trajectories 
for e > 0. 

3.2. Transversality conditions on restriction and lifting. One assumption 
on the restriction 1Z and the lifting C is that they are both smooth maps. 

Furthermore, we assume that the lifting operator C maps some bounded open 
set dom£ C M. d into the basin of attraction U of Cq for e — 0. We will make all 
convergence statements in this section for x € dom C. 

We formulate the transversality conditions on 1Z and C with the help of the 
tangent space Ao(m) to the slow manifold Co in a point uq S Co, which is given as 

Wo(wo) = ker<9/ (uo), (3.6) 

where f E is the right-hand side of the full system ( |2.1[ ), defining the flow M e , in e = 0. 
The definition of the tangent space Ao can be extended to the neighborhood U of the 
slow manifold Co, for example, by using 

Nv(u) := NMu)) (3.7) 

for u £ U that are not in Co- Remember that the stable fiber projection go maps all 
u E U onto the slow manifold Co- The tangent space N e (u) to the perturbed slow 
manifold C E in a point u € C e is a perturbation of Ao(u) of order e. 
Assumption 2 (Transversality of 1Z and C). 

1. The map go° C is a local diffeomorphism between dom£ C K d and the slow 
manifold Cq for every x G dom£. 

Equivalently, the concatenation of the linearizations dgo(C(x)) € M £>x£> and 
dC(x) e M. Dxd has full rank for all x G dom£ C R d . 

2. The map 1Z : U — > M. d , restricted to the slow manifold Co, is a local dif- 
feomorphism between Co and M. d for every u in some relatively open subset 
domftnCo. 

Equivalently, the dimension of the space cflZ(u)J\fo(u) equals d for every u G 
dom 1Z. 
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3. The set dom 7?. n Co contains go(C(domC)) as a subset, and the boundary of 
dom7\!.nCo has a positive distance from the boundary of go(£(dom£)). 
Note that the points [Tl and [2] of Assumption [2] are generically satisfied in a given 
x G R d or uq eW. Point^l] can be visualized meaningfully only if D — d > 2 and is 
illustrated in Figure |3.1[ b) . We define the domains dom £ and dom 1Z such that we 
exclude all non-transversal points. By convention we keep dom£ and dom 1Z such 
that the transversality conditions are uniformly satisfied in dom£ and dom TZ. The 
assumption that dom 1Z (the region where 1Z satisfies transversality assumption [2| 
contains the set go(£(dom £)) guarantees that the map x i-> 7Z(go(£(x))) is locally 
invertible for all x G dom£ and that the linearization is uniformly regular in dom£. 
All points of Assumption [2] and the bound (3.3) on the derivative of the slow flow, 



d\M £ (t\ •) restricted to the slow manifold C £ , can be combined to ensure that the map 

R d D dom£ 3i4 1Z{M £ {t-g e (£{x)))) G R d (3.8) 

is locally invertible for all e G [0, £<j) and t satisfying \t\ < T up /e, where T up i s an 
upper bound determined by the size of dom 1Z. All components of the map (3.8) are 
locally invertible: g e o £ ; dom£ — > C £ by Point [T] of Assumption [2] (transversality of 
£), M £ (t; •) is a diffeomorphism on C e , and 1Z, restricted to Cq (and, hence, to C £ ) is 
also locally invertible due to Point [2] of Assumption [2] For e — the map ( 3 



is 



independent of t. Moreover, the norm of the derivative of the local inverse of the map 



(3.8) is bounded uniformly for sufficiently small e and for x G dom£. 



3.3. Map of exact flow M e into M. d . Next, we give a coordinate system and 
a constructive procedure that maps the flow M E , restricted to the slow manifold C e , 
back to M. d . This kind of map is called a "lifting" of the flow M e on C £ to its cover 
M. d in, e.g., [9], but we do not use this term here to avoid confusion with the lifting 
operation £, used in an equation- free context (cf. for example [23]). For any fixed 
iskip the following map X e : dom£ — > C £ represents (part of) C e as a graph over 
dom£: 

X £ (x) = M £ (t skip ;g e {£(x))). 

This map is locally invertible because go o £ is a local diffeomorphism between dom £ 
and Cq (and, hence, g £ o£ is a diffeomorphism between dom£ and C £ for small e), and 



Af e (i s ki p ; •) is a global diffeomorphism on C £ (see (3.3)). Moreover, if i s ki p < T up /e 
then one can find, for a given u = X £ (x) G C e , a pre- image of any point u G C £ close 
to u by solving 

1Z(X £ {x))=1Z(u) (3.9) 

for x. This follows from Assumption [2] In particular, point [3] of Assumption [2] gives 
the bound on the range of t s ki P for which the linearization of ( |3.9| is regular (the 
trajectory t h-> M £ (t; g £ (£(x))) should not leave domlZ for t G [0,£ s ki p ])- By requiring 



x x, the pre-image x, defined by (3.9) becomes unique. In the coordinate x, mapped 



to the slow manifold C £ by X £ , the flow M £ on C £ satisfies the implicit ODE (see also 



(2.7)) 



^lZ(M £ (t s]iip ;g £ (£(x)))) = ^^(M s (t skip + 6; g s (£(x)))) j^, (3.10) 

whenever e is sufficiently small, 1Z and £ satisfy the conditions listed in Assump- 
tion [2j and for any fixed £ s ki P ■ For different values of t s ki p we get different coordinate 
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representations of the same flow, all related to the representation with i s ki P = via 
the global diffeomorphism M £ (t skip ; •) on C £ , which is a near-identity transformation 
if £ s ki P < 1/e (see (pO|)). 

Let us denote the flow M £ on C s , mapped back into dom£ via X £ , (which is 
generated by the ODE ( |3.10| ) as $*(<5; •) : dom£ -> dom£. If < <5 < 1/e, this flow 
map < &*((5; •) is also defined implicitly by solving the following system for y* : 



TZ(XM) =K(X e (M e (5;x))), 
K(M £ (t skip ;g £ (£(y«)))) =TZ(M e (t skip + 5; g £ (£(x)))), 



that is, 



(3.11) 



and setting <J>*(5;x) := y*. The local invertibility of X £ guarantees that there is a 
solution y* close to x and that the solution y* is unique in the vicinity of x. For larger 
5, one breaks down the flow into smaller time steps such that one can apply the local 
solvability at every step: 



$*{6;x) = $*(8/m;-) m [x] 



(3.12) 



for sufficiently large integer m. This construction of achieves a representation of 
the exact flow M £ restricted to C £ that is globally unique on dom £ for all £ s ki P and 8 
satisfying < 8 < T up /e - t skip . 

3.4. Convergence. We can now compare the result of the approximate flow 
map, y = &(8;x) to the exact solution y* given in (3.111 and compute how the 
difference y — y* depends on x, t s kip 7 8 and s. To highlight where the difference 
between y and y* comes from, we repeat the defining equations of both (implicitly 
given) quantities: 



K(M £ (t sMp ; £(y))) = K(M £ {t skip + S; C(x))), 
TZ(M £ (Ui P ;g e (£(v*))))=K(M e (iskip + 6;g £ {C(x)))). 



(3.13) 



The difference between y and y* follows now from a regular perturbation argument 
comparing solutions of the two equations in (3.131. We rely on (3.4), which guaran- 
tees that the perturbations are small, and the invertibility of the map (3.8), which 
guarantees that the linearization of the left-hand side with respect to y is an order-e 
perturbation of a diffeomorphism. 

Theorem 3.1 (Convergence of approximate flow map). Let K E (0,Kq) be a 
given constant. We assume that the assumptions on time-scale separation (Assump- 
tion^ and transversality (Assumption^ hold for £, M £ andlZ such that 

x^K(M e (t;g e (£(x)))) 

is a local diffeomorphism if \et\ < T up with some T up > that is uniform for all e and 
all x € dom £. 

Then there exist a lower bound t for t sklp , an upper bound Eq for e and a constant 
C > such that y — Q(5;x) and y* = $>*(S;x) are well defined by (3.13), and the 
estimate 



holds for all orders j G {0 
SG [0,T up /e-t skip ]. 



m8;x)-%$.(8;x)\\ < Ccxp(-irt sld p) (3.14) 
,k}, all x G dom£, e G (0,£q)> ^skip G (t Q ,T up /e], and 
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The proof of Theorem 3.1 splits the error &(S; x) — $*(<5; x) using the fiber pro- 



jection g E . The projection of the error onto C e using g e is zero by construction, and 



the error transversal to the manifold decays exponentially due to (3.4). The details 
of the proof are given in Appendix |A"} 

We note that the magnitude of S and t s ki P is only limited by T up /e, where T up 
was defined by the boundaries of domT?. and dom£. According to the error estimate 



(3.14) it would be beneficial to choose i s ki p as large as one can afford computationally. 



If one chooses i s ki P of order l/e p (with p € (0, 1)) then the right-hand side of estimate 



(3.14) is small beyond all orders in s. 

4. Discussion of the general convergence statement and its assump- 
tions. Theorem |3.1| is a local statement with respect to x, claiming convergence only 
in a region dom£ in which the transversality conditions are uniformly satisfied, and 
restricting the times i s ki P and 6 such that the slow flow M e (t; g e {x)) cannot leave the 
slightly larger region dom.ll for the times t = i s ki p and t = i s ki P + S. This is appro- 
priate because in many cases, during continuation or projective integration the maps 
1Z and £ get adapted (for example, for the traffic problem investigated in Section [6j 
£ is varied along the curve of macroscopic equilibria). 

4.1. Comparison to the explicit equation-free approach. The Conver- 



gence Theorem 3.1 implies that a longer healing time i s ki p always reduces the devia- 
tion from the true flow as long as the transversality conditions are uniformly satisfied. 
This is in contrast to the approach, proposed by 23 , where the map $ was explic- 
itly defined as $(6;x) = P £ (S;x) = K(M e (S; £(x))) or <f>(6;x) = P e {8 + i skip ; x) = 
lZ(M E (5+t s ki p ; £{x))) [551I271IS]. Following this approach one would analyse equilibria 
of the slow flow and their stability by studying fixed points of the map 

x^K(M £ (t;£(x))) (4.1) 

for x where 0<t<Cl/eis chosen such that it includes a healing time i s ki P (t > i s ki P )- 
For 1 « ( < 1/e the map 3>(t; ■) : x i— > TZ(M e (t; £{x))) is a perturbation of order 



O(et) < 1 of the map x H> lZ(go(£(x))) . Thus, the explicit approach (4.1) converges 
for e — >• and et — !• only if 7Zog o£ equals the identity on M. d . Often this requirement 
is approximated by 1Z o C = I because g is in general unknown [3TJ [57J |33J |35] . 
Note that there is no e or t dependence in the limiting map 1Z o g a o C, resulting in a 
much more restrictive convergence condition for equation-free analysis based on the 
explicit map <i>(i; x) = 1Z(M E {t] C{x))). In particular, if the map IZo g o £ has robust 
non-trivial dynamics (say, a stable fixed point) then this dynamics will persist under 
small perturbations. Thus, the equation-free approach based on the explicit time-t 
map $(£; lZ(M e (t; C(x))) will generate the dynamics prescribed by IZo g Q o £ 

independent of the properties of the flow on the slow manifold. 

One way to guarantee that lZog o£ = I is to ensure that the lifting operator maps 
onto the slow manifold C e . This has been achieved up to order e through constrained- 
runs corrections to £ |42[ I43j . In our notation the first-order version of this scheme 
would correspond to defining the lifting £ : R d 3 x H> u € R D as the (locally 
unique) u satisfying lZ(u) = x and d / dt(lZ^ (u)) — (zero-derivative principle), 
where 72* is an arbitrary operator satisfying R D = ker TZ © ker "R.™ . Zagaris et al. 
[42], 143] developed general mth-order versions of this scheme. Vandekerckhove et al. 
[39] compared the constrained-runs schemes from (42] 03] to the results of the implicit 



expression (2.9) (called InitMan in 39J ) for various examples, finding (2.9) uniformly 



vastly superior in terms of convergence and performance. Equation (2.9 1 also requires 
only the solution of a d-dimensional, not a Z?-dimensional, system (usually d <C D). 
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4.2. Testing the transversality conditions and choosing the healing time 
and coarse dimension. The conditions listed in Assumption [T] and Assumption [2] 
contain terms that are unknown in practice. For example, the fiber projection go 
and the tangent space Ao to the slow manifold are both inaccessible because in many 
cases one can not vary the time scale separation parameter e. However, observing the 
minimal singular value of the linearization (^^(iskip! x ) = d / dx\TZ{M £ (t s \^ v \ C(x)))] 
with respect to x (a d-dimensional matrix) serves as an indicator: in points where the 
transversality condition is violated, the linearization becomes singular. 

Similarly, the condition of the linearization ^^(iskip! x), cond ^^(^skip! x), is a 
guide how to choose an optimal healing time i s ki p . While the error due to finite time 
scale becomes smaller at a faster rate than cond(92-P e (iskip; x) grows, other errors may 
become dominant when they are amplified by condc^-Pel^skip; x). Especially, when the 
microscopic system is a Monte-Carlo simulation, a trajectory M £ (t;u) is determined 
via ensemble runs, and the accuracy of the evaluation of M £ is only of the order of 
1 / yS where S is the ensemble size. 

The linearization c^-Pe (iskip! x ) a l so helps to discover if one has too many coarse 
variables, that is, if d is too large such that the flow M e restriced to the assumed slow 
manifold C £ is not sufficiently slow (still containing rapidly decaying components). 
Then c^-Pe^skip; x) becomes close to singular, too. Note that any solution found, for 



example, by solving the fixed point equation (2.5 1 is still a correctly identified fixed 
point with correctly identified stability. However, the linearization of (2.5) becomes 
close to singular. 

4.3. Stochastic systems. Barklcy et al. [5] analysed how the equation-free ap- 
proach can be used to analyse moment maps of stochastic systems or high-dimensional 
chaotic systems that converge in a statistical mechanics sense to low-dimensional 
stochastic differential equations (SDEs). They observe that the healing time strongly 
influences the results. Even the inclusion of additional macroscopic variables (in- 
creasing d) drastically changed the results of the equation-free analysis. It is unclear, 



how the implicit scheme (2.3 1 behaves in the situations studied by {2J. While [2] also 
invoke a separation-of-time-scales argument to study approximation quality, their 
setting does not fit into the assumptions underlying Fenichel's theorem, but requires 
weaker notions of convergence (see |13j for a review). An adaptation of the analysis 



in [2] and, possibly, further adaptation of the implicit scheme (2.3), is the missing link 
between the convergence theorem for the idealized situation, given in Section |3j and 
applications of equation-free analysis to stochastic systems. 

5. Traffic Modeling — The Optimal Velocity Model. We now turn to 
the equation-free analysis of a system that illustrates perfectly the framework of the 
general convergence theorem |3.1| We consider N cars driving around a ring road of 
length L. The individual drivers' behavior is assumed to be uniform and deterministic, 
modelled by the optimal velocity model J] of the form 

rx n + x n = V(x n+ i - x n ), 71 = 1,2, ...,N (5.1) 

where x n is the position of car n, r is the inertia of the driver and car, and V is 
an optimal velocity function, prescribing the preferred speed of the driver depending 
on the distance to the car in front (the headway). The ring road implies periodic 
boundary conditions in space 



X n +N =X n +L. 
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(5.2) 




Fig. 5.1. Optimal velocity function V defined in i5M for h = 1.2. The maximal velocity 
r>o(l + tanh(/i)) is obtained for x — > oo. vq acts as a scaling parameter for V. The inflection point 
of the optimal velocity function V is at h. 



In order to do numerical bifurcation analysis the second order ODE (5.1) is rewritten 
as a system of first order ODEs: 



Vn = t~ X [V(x n+1 - x n ) - y n 
Similar to [H [10] we choose the function 



V(Ax n ) — v (ta,nh(Ax n — h) + tanh(/i)) 



(5.3) 



(5.4) 



as the optimal velocity function. In (5.4) «o(l + tanh(/i)) is the maximal velocity, 
Ax n :— x n +i — x n is the headway, and the inflection point h of V determines the 
desired safety distance between cars. The reviews [29l [III [31] put behavioral models 
based on optimal velocity functions into the general context of traffic modeling, and 
discuss possible choices of optimal velocity functions. One conclusion from |31j is 
that the choice of V does not affect the overall bifurcation diagram of a single jam 
qualitatively (some choices of V can give rise to unphysical behavior such as cars 
briefly moving backwards, though). Depending on parameters and initial conditions, 
the system either shows free-flow behaviour, that is, all cars move with the same 
velocity and headway, or it develops traffic jams, which means that there coexist 
regions of uniformly small headways and low speeds, spatially alternating with regions 
of free flow with uniformly large headways and large speeds. We focus on the dynamics 
near the formation of a single jam. After an initial transient, the single traffic jam 
moves along the ring with nearly (due to a finite number of cars) constant shape and 
speed as a travelling wave against the direction of traffic. 

5.1. Direct Simulations. First note, that the uniform flow, starting from initial 
condition 



x n (0) = (n-l)~ 



yn(o) = v 



N 



(5.5) 



is a solution to (5.3), where all cars move with the same velocity y n (t) = V(jj) and 
headway Ax n (t) 



jj. We focus on two types of long-time behavior, the uniform 
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(a) h = 1.2, d = 0.87 (b) h = 1.2, v = 0.91 

Fig. 5.2. Comparison of the two traffic flow regimes. The initial condition (blue) is compared 
with the final state (red = T = 5 X 10 4 , green= T — 500). (a) Free flow regime, (b) Traffic jam 
regime. Note the different scales of Ai„ on the vertical axis. 




(a) h = 1.2, D = 0.87 (b) h = 1.2, v = 0.91 

Fig. 5.3. Time evolution of the macroscopic variable a for the same parameters as in Figure 
\5.S\ (a) shows the decay to the stable free flow. (6) Using the same initial condition as in (a), the 
system converges to a stable traffic jam. The inset in (b) shows the difference to the average value 
ct* . One expects small oscillations of a in time due to the finite number of cars. However, these 
small oscillations are below the tolerance of the ODE solver. 



flow and travelling wave solutions. To give a qualitative picture of these, we run two 



simulations, initializing the system (5.3) with initial conditions close to the uniform 
flow 



... L ( 2tt 

Xn(0) = {n- 1)— + /ism I — n 



Vn(0) = V 



(5.6) 



where /x is the strength of the periodic perturbation. For all simulations, we use 
N = L = 60. The simulations were run for a time T = 5 • 10 4 using Matlab's 
ode45-solver jT8j with absolute and relative tolerance 10 -8 . All parameters for the 
simulation can be found in table [B] in Appendix [Bj In order to compare simulation 
results, h = 1.2 has been fixed. Furthermore, the velocity parameters vq — 0.87 and 
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Vq = 0.91 have been chosen to illustrate the uniform flow and the traffic jam regime, 
respectively. In Figure |5.2[ the headway is shown as a function of car number. It 
can be seen, that the initial perturbation either decays to the uniform flow or the 
trajectory converges to a travelling wave solution, depending on the choice of vq. 

We choose the standard deviation a for the headway as the macroscopic measure 
(called x in Section [2] and Section |3| describing the traffic flow 



\ 



N - 



1 N 



(Ax)) , where Ax„ = x n +\ 



(5.7) 



Here, (Ax) = 1 is the mean of all headways. 



and the decay of a to the free flow is shown in Figure 5.3(a) 



The free flow corresponds to a = 
If vq is chosen in the 



traffic jam regime, a increases until it settles to an equilibrium, where a travelling 
wave of fixed shape is observed. It can be seen in the inset of Figure [5.3(b) that the 
macroscopic variable oscillates even in its steady state. These small-scale oscillations 
are expected due to the finite number of cars, because cars arrive at the rear and leave 
at the front of the jam at periodic intervals. However, the oscillation amplitude is 
orders of magnitude smaller than the macroscopic dynamics such that the oscillations 
are obscured by discretization effects of the ODE solver (which shows sub-tolerance 
oscillations even for systems with stable equilibria). The next section presents an 
equation- free bifurcation analysis for the jam formation on the macroscopic level. 

6. Equation-Free Bifurcation Analysis. We choose a one-dimensional macro- 
scopic description, that is, the standard deviation a is the only macroscopic variable, 
even though we expect the overall slow manifold to consist of several patches, some of 
which are higher dimensional. The other patches correspond to states with coexisting 
traffic jams, which are meta-stable (transient, but existing for extremely long time). 
Coexisting jams move with a speed difference of order exp(—C/N p ) for some C > 
and p > [31] and take exponentially long (in N) to dissolve or merge. 

The change of the chosen macroscopic variable er is studied with respect to sys- 
tem parameters. According to the equation-free approach presented in Section [2] the 
macroscopic ODE has the implicit form 



^-K{M(t skip ,£{a)) = ^K(M(t skip + 6,£(a)) 



8=0 



(6.1) 



where the derivative on the right-hand side is approximated by the finite-difference 
quotient with finite 5 



F(a) 



KM(t skip + 5, C(a)) - KM(t skip , C(a)) 



(6.2) 



and t s kip is the healing time. 

As explained in Section [2] and |3j the equation-free setup avoids an analytical 
derivation of a macroscopic ODE but uses (6.1 ) where (6.2 ) is evaluated by simulation 
bursts of length i s kip + ^- Note, that M and C depend on the system parameters h and 
vq, which are not expressly included in (6.1) and (6.2 1. We also drop the subscript 



e of M because it does only enter our system indirectly through the number of cars 
N. In order to find trajectories or equilibria of ( 6.1 1— ( 6.2 ) , it is necessary to define a 
lifting operator C and a restriction operator 1Z. In our case, the restriction operator 

16 



7Z is given by the definition of the macroscopic measure in eq. (5.7), i.e., 



K(u) 



\ 



N 



1 N 



(Ax)Y 



(6.3) 



Our lifting operator constructs initial conditions with the help of a reference state 
u = (x,y) £ K 2Ar , obtained during a microscopic simulation. We have to guarantee 
that the lifting C initializes the system into the vicinity of the solution of interest. 

This requirement was listed in Section [3] as C having to map into the attracting 
neighborhood U of the slow manifold. 

The following description concentrates on a single-pulse traffic jam. The compo- 
nents of the reference state u are the positions (x n )^ =1 and the velocities {y n )n=i of 
the cars, respectively (cf. ( |5.3[ )). Let us denote the macroscopic state corresponding 
to the restriction of u by a = lZ(u). The lifting operator is then defined as 



C p (u, a) =u = (x, y) = (x ne w, 2/ncw) eK x 



*£new. 1 



0. 



N , where 
} ^(Ax-(Ax))+{Ax), 

n-l 

^ Asncw,! n = 2,...,N, 



(6.4) 



V(Ax ncw , n ) n = l,...,N. 



V is the optimal velocity function (5.4), (•) refers to the average of a quantity, and Ax 
are the headways of the reference state (Ax n = x n +i — x n ). In (6.4) we compute the 



positions x £ K first and then initialize the velocities y £ K by using the optimal 
velocity function for these positions. The positions are initialized such that x\ = 0, 



resulting in a unique mapping from headways to positions. The definition (6.4) of £ p 
contains an artificial parameter p, which we keep equal to unity throughout, except 
for Figure [6~2] in Section [6TT] A parameter value of p ^ 1 introduces a systematic bias 
into our lifting such that we can vary p gradually to investigate how our results depend 
on our choice of lifting. For p ^ 1, the lifting £ p violates the common assumption of 
equation-free computations, where the identity 1Z o C — I is claimed to be necessary 
[531 12"61 |2"T1 |2"T1 13"6] . An application of C and 1Z without any time evolution in between, 
yields TL{C p {a)) — p- a. 

In the following, we use an equation-free pseudo-arclength continuation scheme 



to compute bifurcation diagrams for the fixed point of (6.1 )-(6.2), that is, we track a 
root curve (branch) of 



F(a,v ) = Q 



(6.5) 



in the (a, v )-plane. The influence of speed limits on traffic jam formation motivates 



the choice of the velocity parameter vq as a bifurcation parameter. In (6.5) we include 
the bifurcation parameter vq explicitly as an argument of F. The pseudo-arclength 
continuation contains two steps. The first step is a predictor step, where we use a 
secant predictor assuming that we know two points on the branch already. Let (er°, Vq) 
and (cr 1 , - ;^) be those two points. We define the secant direction by 



(- 1 



v° ). 



(6.6) 
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The prediction (<7,t»o) for the next point on the branch is then determined by the 
secant predictor 



(a, Ho) = (o- 1 ,^) + s "pj|' ( 67 ) 

where we keep the stepsize of the predictor uniformly at s = 10~ 3 . The prediction 
is not exactly on the branch and must be corrected in the following corrector step, 



which is chosen to be perpendicular to the predictor direction (6.6 1. The corrector 
step solves the system 

,/ (< "*° , = <"> 

w {a) (a - a) + w {vo) (v -v ) = 0, 



where uA CT ) and w^" ' are the components of w, respectively. System (6.8 1 can be 
solved with respect to a and vq by Newton's method using 

(a k+ \v% + y = (a k ,v$) + vF^F{a k A), (6.9) 



where Fx is the Jacobian of the left-hand side of (6.8), given by 



and v is a relaxation parameter adjusting the length of a Newton step. For all compu- 
tations we used a full-Newton step, that is, v = 1. If the information on the Jacobian 
of the system is poor, for example, in noisy or stochastic systems, it might be useful to 
use a damped Newton method (v < 1). The iteration is initialized with the predictor 



(6.7) 



(a , v° ) = (*,%). (6.11) 
During the iteration the function F has to be evaluated always by lifting, running 



simulations of the microscopic system and restricting, according to its definition (6.2 ) 
with t skip = 300 and S = 2000. 

The Jacobian Fx is approximated via finite differences. Since w^ a ' and if/ 1 " ) are 
known from the predictor step, we only have to determine F a and F Vo . We evaluated 
F at the points 

(a, w ), (er + Aer, v ), {a, v + Av ) (6-12) 
and computed the one-sided derivatives 



F(a + Aa,v )-F(a,v ) „ F(a,v + Av ) - F(a,v Q ) 
F a = ^ , F V0 = _ (6.13) 



We started the one-parameter continuation of the traffic jam in direction of de- 
creasing vq from two profiles obtained by direct simulations at v — 0.91 and v = 0.9, 
respectively. The resulting bifurcation diagram is shown in Figure |6.1| 

The traffic jam, i.e., travelling wave, is stable for large values of vq- When fol- 
lowing the branch, a saddle-node bifurcation is detected at (1Z(M(t s ki p , Ca))* , Vq) « 
(0.88,0.125), where the traffic jam changes stability. A further decrease of vo at 
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0.88 0.881 0.882 0.883 0.884 0.885 0.886 0.887 0.888 0.889 0.89 

v o 

Fig. 6.1. Bifurcation diagram obtained by equation-free pseudo-arclength continuation for h = 
1.2. The traffic jam profiles are shown for selected points in the bifurcation diagram, marked with 
black circles. Note the change in scale on the vertical axis on the profiles for a better visibility 
(horizontal axis shows the car number n, vertical axis show headways) . A fold point has been detected 
at (vq, 1i(M (t s kip, "£"■))) ~ (0.88,0.125), where a change in stability is observed. The blue dots mark 
stable states, while the red dots mark unstable states. It is due t o th e equa tion- free continuation, 
that unstable branches can be observed. For the lifting we use | |5.6| and | |6.4[ l for continuation 
of the uniform flow and the travelling wave solution, respectively. Additionally, the black crosses 
mark a backward trajectory computed by using j6,14| . Starting from the stable branch, the backward 
integration converges to the unstable branch (big cross). The black dot is the base point used for an 
error estimate in i6.20\ . 



that point would make the traffic jam dissolve. But due to the equation-free pseudo- 
arclength continuation of the continuous branch, it is possible to follow the branch 
around the fold point and continue the unstable branch for increasing vq. The traffic 
jam stays unstable until it reaches the uniform flow at a = at a Hopf bifurcation 
point (cf. Section 6.3 and (6.26)). The microscopic states corresponding to selected 
points along the branch are shown as insets in Figure [6TT) The shape has sharp layers 



and a flat plateau on the stable branch, and becomes harmonic close to the equilib- 
rium value a = 0. Additionally, the time steps of a backward integration are shown 
for vq = 0.884, showing the heteroclinic connection between stable and unstable jam. 
The trajectory starts for to — at the stable branch and the Euler scheme (2.8) is 
used for computing the backward trajectory, that is, 



HM(t skip ,C(a j+ i)) = nM(t skip ,C(aj)) + F(a 3 )At, 



(6.14) 



where it, is the solution at tj — jAt, and At = —5000 is chosen. For the compu- 
tation of F(a) the parameters from table |B| are chosen in eq. (6.2). The backward 
integration converges to the unstable branch. We will use this backward integra- 
tion in combiniation with a forward simulation to obtain an error estimate for the 



equation- free approximation in Section 6.2 and (6.20) 
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Fig. 6.2. (a) Bifurcation diagram obtained from running the implicit equation-free continua- 
tion scheme described in Section^ for p = 0.95,0.97,0.99,1.0,1.02,1.05. The black dots show the 
results from a direct downsweep of the stable branch. Depending on the value of p, the results differ 
significantly from the direct simulation data, which is usually used as an argument to find a 'good' 
lifting operator, (b) The healed version of the bifurcation diagram, i.e., a"h ca iod = TZM (i s ki p , C p a) 
for different p all collapse to the same curve, fitting the direct numerical results perfectly. See also 
the main text in Section^ (c) Analysis of the lifting error. The blue data points show the distance 
between the equation-free solution vo(cr) and the restriction of the simulation data using eq. (6.16) 
as a measure for the error. The results from a direct simulation of the stable branch are used as 
a reference curve (cf. Figure \6.2(a)\ . For the 'normal' equation-free data, i.e., using the unhealed 
macroscopic quantities, it is observed, that the error is minimal at p = 1.005, corresponding to a 
'good' lifting operator. The g reen da ta points show the behaviour of the healed version of the bi- 
furcation branches (cf. Figure \6 '.2(b)f . No significant error can be observed when using the healed 
data. 



6.1. The influence of the choice of lifting operator. Figure 6.2 shows how 
the results depend on the artificial parameter p, which we introduced into the lifting 
operator C p . In both panels, the same bifurcation diagram is shown for several values 
of p and compared to the restrictions of the stable fixed points of direct long-time 
simulations (T = 3 • 10 5 , black dots). The case where the usual equation- free identity 
1Z o C p — I is fulfilled corresponds to p — 1. We observe that the pre-images a of the 
equilibria under the combination of lifting operator and healing M(t s ki p ; jC p (-)) depend 
significantly onp (panel (a) of Figure 6.2). Therefore, we compare Figure 6.2(a) with 



the corresponding Figure 6.2(b) for the healed macroscopic quantity 



Chealcd = ~R(M (tg^ip, £ p a)) 



(6.15) 



for each macroscopic equilibrium a along the branch of the bifurcation diagram. Ac- 
cording to Section [3] the map lZ(M(t s ki P , C p cr)) is a local diffeomorphism from R d 
into M. d with d = 1. Plotting the bifurcation diagram in the (vq, CTheaied)-plane in 



Figure 6.2(b) we obtain a solution branch that is independent of the choice of the 



lifting operator, as expected. 

For a more detailed analysis of the error, we compute the L2 norm between the 
interpolated data sets vq(<j) (expressing the parameter as a function of the equilibrium 
location near the fold) for the direct simulation data and the data for the stable 
branch of the equation-free bifurcation diagram. For interpolation, Matlab's interpl 
function [T5] with the ' spline ' option is used. We use the error measure 



11/ -.911 



[f (a) - g(a)} 2 da, 



(6.16) 



to analyze the deviation between the restriction of the direct simulation data and 
equation- free continuation data. Here, / and g are the interpolated data sets vq(ct) 
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(a) (b) (c) 

Fig. 6.3. (a) Bifurcation diagrams for h = 1.12 and t s ki p = 300 (blue curve) and i s ki p = 
2000,4000 (green curve, red curve), respectively. The difference between the curves is very small, 
(b) Comparison of the Jacobians for the different values of t s ki p along the curve. A significant 
change in behaviour is observed for different values of t s ki p . (c) The dependence of the result on 
<skip- The origin on the vertical axis is arbitrary since the correct solution of the full system is 
unknown on the unstable part of the branch. We set the y-origin at the result for t s ki p = 2000 (red 
dot). 



for the simulated data and the equation-free data, respectively, in the range of a 
between a — 0.125 and b = 0.25. The unstable branches can not be compared 
with direct integration of the system. The deviation E using eq. (6.16)) with lifting 



parameter p is shown in Figure 6.2(c) The blue data points correspond to the distance 
between the restriction of the simulation data and the equation-free solutions (that 
is, the pre- images of the equation- free microscopic solutions under M(i s ki p ; £ p (-)) in 
the domain of C p ) . The distance is small for p values close to 1, where the usual 
identity 1Z o £ = I is fulfilled. However, the distance for <7h ea ied (green data) is 
uniformly small, independent of the choice of p. Therefore, healed quantities should 
be used when comparing equation-free results to restrictions if the direct simulation 



data. The uniformly small errors in Figure 6.2(c) (in green) suggest that with implicit 
time steppers the results are not sensitive to the choice of the lifting operator. This 
is in contrast to most equation- free applications 23, 5, 21 , which used explicit time 
steppers of the form $(5; x) — TZ(M (S; £(#))). 

6.2. The influence of the healing time i s ki P - In this section, we investigate 
the influence of i s ki P on the equation- free results, e.g., bifurcation diagrams and sta- 
bility analysis. First, we show that the bifurcation diagrams are rather insensitive to 
the choice of i s ki P , while the information of the Jacobian depends more noticably on 
the value of i s kip- 

The bifurcation diagrams obtained for h = 1.12 and £ s ki P = 300,2000,4000 are 



shown in Figure 6.3 respectively. In Figure 6.3 (left panel), it can be observed, that 



the bifurcation diagrams are similar for all choices of i s kip, i.e., they show the same 
qualitative features. Although the bifurcation diagrams are quantitatively close to 
each other, the information about the derivatives, i.e., the Jacobian dF/da, differs 
significantly. For i s ki p = 300, the fold points does not correspond to a sign change 
of the Jacobian dF/da (cf. Figure 6.3 1. The approximate Jacobian is monotone 
increasing, approaching zero from below for a — > 0. In particular, the Jacobian does 
not change its sign near a — 0.07, the upper fold in Figure |6~3^ a). Using a longer 
^skip = 2000 the fold points are detected by Jacobian sign changes. Using an even 
higher t s ki p = 4000, it can be observed, that the Jacobian changes sign only once, 
yielding an inconsistent behaviour, again. Therefore, the Jacobian depends on £ s ki p 
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Fig. 6.4. Error by t s ki p and At = —28 when using the forward-backward integration as a com- 
parison. The color is a logarithmic coding (to base 10) for the error according to | |6,20| . Projections 
into the error-parameter plane can be seen in Figure \6. 5\ 
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Fig. 6.5. (a) Error by ( |6,20| l for varying At = —28. Data is plotted on logio axes, yielding a 
power-law dependence of the error on (At) 2 . The black line indicates a power law with exponent 2 
as reference. This result suggests using At and 8 larger than 300 to be in the region with quadratic 
error estimate. For the continuation presented in Section^ we used 8 = 2000. (6) Same error- 
projected to igkip ■ The error increases with At = —28, but is roughly constant over the shown range 
of *skip ■ 



substantially. In order to investigate the error caused by t s ^ v systematically, we use 
the bifurcation diagram obtained for £ s ki P = 2000 (which has the most accurate sign- 
change of the Jacobian compared to the actual bifurcation diagram) as a reference 
and show the difference between bifurcation diagrams for other values of t s ki P and 
this reference. For different values of t s ki P the distance to the reference bifurcation 
diagram is measured by using 

\\f-9\\ = [ \f(a)-g(a)\da, (6.17) 

J a 

where / is the reference diagram for i s ki p = 2000 and g is the diagram for other 
values of i s ki P - We exploit that the diagram is a graph over a near the folds and use 
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The results are shown in Figure 
min) of the vertical axis gives a good estimate for the order 
of magnitude of the error for the bifurcation diagram. 

Another way of determining the error caused by the finiteness of i s ki P is to use 



a = 0.01 and b = 0.28 as interval boundaries for a 
The range (max 



6.3(c) 



the equation-free backward integration scheme (6.14). We use a point a(t) on the 
macroscopic backward trajectory (black dot on the backward trajectory in Figure 
6.1) as a starting point for our computation. From this point, we compute a healed 
value ctq by 

nM(t sUp ,La(t))=:a^ (6.18) 

This point is used as the reference solution. Then we apply a backward Euler step 

At) =: <7i (using a Newton method with 
From cri we run the microscopic simulation for i s ki P — At, yielding 



(6.14) to find the implicitly defined a(t 
7 ) 



tolerance 10 
the value 



ftM(i skip -At,£ai) 



(6.19) 



(note, that At was chosen to be negative in the backward integration scheme). The 
values of <7q and Uq should be the identical except for two possible sources of error, 
the discretization error of the Euler step of size At and the incomplete convergence 
of the trajectories to the slow manifold during the healing time t s ki P - Therefore, their 
difference can be used as an error estimate: 



error 



(6.20) 

and At = — 2<5 is shown on 



plicit Euler 
gives a 



6.4 



In Figures 6.4 and 6.5 the error for different values ofi s ki P 
a logarithmic axis. For the choice At — —25 the one-step error for the im 
results in an estimate of order (At) 2 . The two-parameter plot in Figure 
qualitative overview about the overall behaviour of the error. Figure [6~5| shows cross- 
sections of this plane. It can be seen, that the error increases like a power law with 
(At) 2 (note the logarithmic axes in Figure 6.5 '&)). For small values At = —26 the 



inaccuracy in the finite-difference approximation (6.2) becomes dominant, yielding a 
deviation from the expected (At) 2 behavior. With respect to i s ki p the error stays in 
the same order of magnitude in the investigated region, such that we conclude that 
the error is rather insensitive to a change of i s ki P over the range we investigated. For 
continuation in Sect ion |6| w e used 6 — 2000 to be in the regime for the quadratic error 
estimate (cf. Figure 6.5(a)). 

6.3. Continuation of the fold in two parameters. A two-parameter scan, 
showing one-parameter bifurcation diagrams in the velocity parameter vq for different 
values of the safety distance h, is presented in Figure 6.6 'a). The curve of folds as a 
result of two-parameter continuation in Figure 6.6 a) shows how the fold merges with 



another saddle-node point in a cusp. The system of equations for continuation of the 
fold is 1251 



with the Jacobian 



F(a,v ,h) = 
F a (a,v ,h) = Q 
) M (a-a) + w (vo) (v -v ) + (h-h) = 



Fx = 



Fa 


r v 


F h 


7(7 


TP 

r v a 


Fha 




w {vo) 


w (h) 




23 





(6.21) 



(6.22) 



h = 1.08 to 1.25 




(a) (b) 

Fig. 6.6. (a) Bifurcation diagrams for h £ [1.08, 1.25], were h increases from the left curve 
to the right curve. (6) Two parameter continuation of the fold point in the parameters vo and h. 
The blue crosses mark the points determined from the bifurcation diagrams of the one parameter 
continuation and the red circles are the results from a two parameter continuation. The black lines 
show analytical results for the Hopf bifurcation at a = (cf. | |6,26| ). The numerical results for 
the Hopf bifurcation (zeros from bifurcation diagrams in the left panel) are denoted as the green 
dots. They are in perfect agreement with the analytical results. Note that in the parameter plane 
projection, the difference between the Hopf and the fold point is barely visible and the first Hopf 
curve and the numerical data is obscured by the numerical data for the fold continuation. 



Since derivatives of second order are needed, we apply an approximation of second- 
order accuracy for the derivatives, i.e., centered differences for the parameter deriva- 
tives in vq and h and one-sided second order schemes for derivatives in a. We use the 
one-sided second-order approximation for F a , because a is non-negative by definition. 
Details for the numerical evaluation of the derivatives can be found in Appendix [C| 

During the two-parameter continuation the Newton iteration used full Newton 
steps [y = 1 in (6.9)). Panel (b) of Figure 6.6 shows the results. They are in perfect 
agreement with the data obtained by a one-parameter continuation. For comparison 
we have included the Hopf bifurcation point of the full microscopic system at a = 0. 
The Hopf bifurcation is a pitchfork bifurcation at the macroscopic level. However, 
since the standard deviation as macroscopic measure is non-negative by definition it 
shows only the non-negative branches. The analytic expression for the Hopf bifurca- 
tion parameter can be found by linearizing system ( |5.3[ ) around the uniform flow and 
using the ansatz (x n (t) , y n (t)) = (x n (0) exp(icjt), y n (0) exp(iuit)). This results in the 
system 



Vn 



iuy n = t 



T' ( — ) {x n+1 - x n ) - y n 



(6.23) 
(6.24) 



where u> is the frequency and V'(jj) the first derivative of the optimal velocity function 
at equilibrium. Eliminating x n and using the periodic boundary conditions results in 




(6.25) 



This implicitly defines vo as a function of h (through V) and can be solved for our 
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specific choice of V (see (5.4)) to yield 



1 - cos(2vrj/iV) 
rsin 2 (2irj/N) (l - tanh 2 (h - £)) 



v ° = _. ; __2 , n __- /^n /i ,.__ 1 _2 77 TTY' ( 6 - 26 ) 



where j = {1,2, . . . AT — 1}. The Hopf curves for the first four spatial frequencies 
(j = 1,2,3,4) are shown in Figure [6.6(b) The analytical results for the first Hopf 



curve are in perfect agreement with the numerical data. Note that the curves for the 
Hopf bifurcation point and the fold point are close to each other in the parameter 
plane shown in Figure [6~6) (b) . 

7. Conclusion and Outlook. In this paper we derived an implicit method for 
equation-free analysis and proved its convergence for slow-fast systems with transver- 
sally stable slow manifolds. We gave a demonstration by performing an equation-free 
bifurcation analysis on a one-dimensional macroscopic description emerging from a mi- 
croscopic traffic model based on a deterministic optimal velocity model for individual 
drivers. We demonstrated that the obtained bifurcation diagrams are independent of 
the lifting operator and the healing time in a suitable region. The bifurcation diagram 
shows a saddle-node bifurcation which is continued in a two-parameter equation-free 
pseudo-arclength continuation. 

The proof of convergence for the implicit coarse-level time-stepper asumes that 
the slow manifold is transversally stable. The review [13] lists in which sense a fast 
high-dimensional chaotic or stochastic system converging in the mean can be viewed as 
a slow-fast system converging to its slow manifold. In practical applications the result 
from Section |3.1| may be used as a plausibility check: the equation- free methodology 
of Kcvrckidis et al appeals to the notions of singular perturbation theory (cf. the 
illustrative example in |24j ) . For any particular system under study, one can check if 
this intuition is indeed justified by testing if the results for the implicit time stepper 



given by ( 2.3 ) are indeed independent of the lifting C and the healing time t s k ip if one 
varies both gradually. This appears not to be the case, for example, for the moment 
maps studied in [2J, as demonstrated in the paper [2] itself. 

For the traffic problems studied in our paper, one long-standing problem is the 
motion of several phantom jams, i.e., multi-pulse solutions, relative to each other. 
For a large number of cars (including the N — 60 cars we used) this motion is 
very slow and therefore near impossible to observe in direct numerical simulations (a 
phenomenon that is called meta-stability) . An open question is whether one can derive 
a computable criterion that predicts for a given configuration of several jams and given 
driver parameters, which of those will collapse or merge and when. This criterion 
might be based on the shape of the traveling wave. One particularly appealing feature 
of equation-free analysis is that one can continue macroscopic equilibria in N, the 
number of cars, using the microscopic model. 

Models closer to situations of practical interest, say with more realistic opti- 
mal velocity functions, randomly assigned driver behavior parameters, an element 
of randomness in the driver behavior or multiple lanes as discussed in the literature 
[T71I2II1I2I] are also amenable to equation- free analysis. This should provide additional 
information to help match parameters of macroscopic models to microscopic driver 
and road parameters. 
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Appendix A. Proof of Convergence Theorem |3.1[ For the proof of Theo- 
rem 



3.1 we have to analyze the two equations (for y and respectively) 



K(M £ {tsM P ; C(y))) = K(M £ (t skip + 6;£(x))), (A.l) 
■R.{M e (t S K P ;ge{£{y*)))) = TZ(M e (t skip + 6;g e (C{x)))). (A.2) 

In both equations x £ R d enters as a parameter. For ( |A.2[ ) we have established 
already in Section [3] that there exists a solution y*, and that it is locally unique. The 
equations ( 3.11[ ) and ( |3.12 ) gave a procedure to pick in a globally unique way by 



starting with = x for 8 — and then extending the solution for varying 5 until 
one reaches the desired value of 5. This procedure achieves unique solvability for 
for all x £ dom£ and for i skip £ [toiT up /e) and S > satisfying t skip + 5 < T up /e. 
For equation (TAT]) we have to prove the existence of a solution y, and prove that it 



is close to (including all derivatives with respect to x up to order k). 

In order to do this, we need to make the consequences of Fenichel's Theorem more 
explicit. The Fenichel result (3.3) implies that the map (U is the neighborhood of Co, 



which also contains C e ) 

T e : K x U 3 (t, u) i-> M e (r/e; g e {u)) 

is well defined and k times differentiablc for all u £ U and all r £ K. This map T e is 
the flow map when restricted to the slow manifold C £ , and projects all points in the 
neighborhood of the slow manifold C e along the stable fibers using g £ . Note that r is 
the time on the slow time scale as we divide by e in the evaluation of the map. The 
derivatives of T e with respect to its second argument u are uniformly bounded for all 
u £ U as long as as r £ [0, Tup]: 

\\diT £ (r;-)\\<C, (j =0...,fc, and r e [0,T up ]). (A.3) 

Correspondigly, the map 

A e {t; ■) = K{T e {r; £(•)) : dom£ 3 x ^ 1Z(M £ (t/s; g e (C(x)))) (A.4) 

is well defined for all r £ K and locally invertible for all e £ [0,£o) an d r satisfying 
|t| < Tup. Note that the range of admissible e includes e = because the limit of 



the right-hand side of (A.4 1 for e — is well defined as the solution of a differential- 
algebraic equation on Cq on the slow time scale. The norms of the derivatives of A e 
and its (locally unique) inverse can be bounded by a uniform constant C independent 
of e £ [0, eo) and r as long as |r| < T up : 

\\diA e (r;-)\\<C, \\diA-\r--)\\<C. (A.5) 

Similarly, the motion transversal to the slow manifold C e consists of a fast decay and 
a slow tracking of the dynamics on C e . Let K < Kq be a given contraction rate, and 



choose the upper bound eo such that the contraction property (3.4) of the stable fiber 
projection g E holds for all e < £o an d all u £lA. Then we can express the transversal 
component of the flow starting from an arbitrary u ElA and t > in the form 



M £ (t;u) = M e (t; g e (u)) + exp 
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(A.6) 



(this defines M^~). In the right-hand side of ( |A.6[ ) the map is k times differentiable 
with respect to its argument u for all e G [0, Eq) (including s = 0), and the norms 
of M £ (t;u) and its partial derivatives M^(t;u) are uniformly bounded for all t G 

[0, oo), e G [0, So) and u &U: 

\\M^(t;x)\\ < C, \\4M^(t;x)\\ < C. (A.7) 

The pre-factor exp(— Kt) can also be extracted if the smooth restriction map TZ is 
applied to both terms on the left-hand side of (A.6), and if we insert £(x) for u: 

K(M £ (t;C(x)))-K(M £ (t;g £ (£(x)))) = 



9R(M e {t;g e (C{x)))+p[M s {t;C{x)) - M £ {t; g £ {C(x)))}) dp 

x [M £ (t;C(x))-M £ (t;g £ (C(x)))] 
dlli^eiet; C{x)) + pexp(~Kt)M^-{t; £(x))) dpexp{-Kt)M^-(t; C{x)) 



(A.8) 
(A.9) 



(A.10) 



We applied the mean- value theorem to equate (A.8 1 and (A.9). To get to the right- 
hand side of (A.10), we inserted the representation (A.6) and used the definition of 
the map T £ . This right-hand side in (A.10) has the form 



right-hand side of (A.10) = cxp(— Kt)r £ (et, t; x), 



(A.ll) 



where the first argument of r £ refers to the time dependence of A e in the argument 
of dlZ. Note that we have introduced the slow time scale as an additional argument 
into r e . We will consider r £ (r, i; x) for arbitrary r G [0, T up ] and t G [0, oo) below, and 
later insert r = et as a particular case. The map r e (r, t; x) is k times continuously 
differentiable with respect to x. The norm of r £ and the norm of its derivatives with 
respect to x are uniformly bounded for x G dom£, e G [0,£o) and r G [0,T up ] and 
t G [0, oo) because all of its ingredients have bounded derivatives (listed in (A. 3), 
dA^7|): 

\\r £ (r,t;x)\\ < C, ||^r e (r,*;a:)|| < C (A.12) 

, k}. Let us define the times corresponding to t s ki P and 5 on the slow 



for j G {1, . . . 

time scale as: 



7"skii 



et 



skip i 



eS. 



(A.13) 



If Tskip and T s ki P + A are in [0, T up ] then the solution of the exact flow satisfies 
(using the locally invertible map A £ defined in ( |A.4[ )) 



A e (r skip ;y*) = A e (r skip + A;x). 
Using r £ and A £: equation ( |A.1[ ) can be re- written as 

A £ {T skip ;y) + sir £ (T s]iip ,t 1 ;y) = A e (r skip + A;x) + s 2 r e (T skip + A,t 2 ;x), 
where 

si = cxp(-X< skip ), s 2 = exp(-if (t skip + 5)), 

tl — ^skip > 0, t 2 = <skip + S > 0. 
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(A.14) 
(A.15) 

(A.16) 



We will first consider solvability of (A. 15 1 with respect to y for general s% and s 2 close 
to 0, and t\, t 2 £ [0, oo). This solution y will depend on the parameters s\, s 2 , t\ 
and t% (among others). Whenever we subsequently insert the particular values from 



(A.13I and (A. 16) for r s ki P , A, si, s 2 , ti and t 2 , the solution y of (A. 15) becomes 
also a solution of ( A.l\ . For each of the terms, A e , A~ x and r E , we have uniform 
upper bounds ((A. 5 1 and( A.12| ) for their norms and all derivatives up to order k for 
the entire range of arguments: x, y £ dom£, r s ki p € [0, T up ], T sk j p + A £ [0,T up ], 
t\, t 2 £ [0, oo), e £ [0, eq) (where £q is determined by the choice of decay rate K as 
given by Fenichel's Theorem). Thus, we can use (A. 14) and (A. 15 1 to establish the 
existence of y, and its distance to y* using the Implicit Function Theorem at the point 
si = s 2 = 0. 



The exact solution y* is a uniformly regular solution of (A. 15) for s% = s 2 = 0, 
all x £ do m£, g £ [0,e ), 7"ski P £ [0,T up ] and A £ [-r skip , T up - r skip ]. Thus, for small 
Si and S2, (A. 15) has a locally unique solution y £ dom£ which depends smoothly on 



all parameters (we write y(x, Si,s 2 ) to emphasize the dependence on (si,s 2 ) £ K 2 ) 
such that 



\\d{y{x,si,s 2 ) - (^(a^Hoo < C||(si, s 2 )||oo 

(j £ {1, ...,&}) for some constant C and all si, s 2 £ (—p,p) for some p > 0. 
Consequently, if we choose to such that exp(— Kto) < p and decrease £q such that 
*o < T up /e then we have for all e e (0,e ), *ski P € [io,T up /e], J € [0,T up /e - i skip ] 
and a; £ dom£ that 



d{y(x,exp(-Kt skip ),exp(-K(t skip + 5))) - 9 J y*(x) 

< C||(exp(-^ skip ),exp(-^(i skip + ( 5)))|| ! 

< Ccxp(-i ; i:f skip ) 



< 



for all j £ {!,..., k}. This establishes the convergence claim of Theorem 3.1 since y 



is the solution of ( |A.l[ ) if si = exp(-KT skip /£) = exp(- Kt skip ) , s 2 = exp(-i4T(T skip + 
A)/e) = cxp(-fsr(t skip + 5)), ti = t skipi t 2 = t skip + 6, r skip = Et skip and A = eS. 

Appendix B. Parameters. The parameters used for the simulations are listed 
in Table El 



Parameter 


Value /Range 




1.7 


L 


60 


N 


60 


A* 


0.1 


s 


0.001 


6 


2000 


At 


-5000 


^skip 


300 




0.8 ... 1.0 


h* 


1.0 ... 1.7 



Table B.l 

Parameters for numerical studies. The quantities marked with asterisk (*) are bifurcation parame- 
ters, where the used range is noted 
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Appendix C. Finite Differences. For the scheme (6.21 ), F is evaluated at the 
17 points 



1 : (a,v ,h), 

2 : (<r + Aa,v ,h), 
4 : (a + 3Aa,v ,h), 
6 : (<r,vo - Av ,h), 

8 : (cr + Act, w — Av , h), 
10 : (cr + 2Aa,v Q - Av ,h), 
12 : (a,v ,h- Ah), 
14 : (cr + Acr, w ,/i - Ah), 
16 : (cr + 2Acr, w , h -Ah), 



3 
5 
7 
9 

11 
13 
15 
17 



(cr + 2Acr, Uq, h), 

(cr + 4 Acr, v ,h), 

(cr, w +Au ,/i), 

(cr + Acr, v + Av , h), 

(a + 2Aa, v + Av a ,h), 

(cr, Vq, h + Ah), 

(a + Aa, v a ,h + Ah), 

(cr + 2Acr, Vq, h + Ah). 



(C.l) 



where Acr = Avo — Ah — 0.001 are offsets for the approximation. One can use 
the following second-order accuracy scheme to compute the derivatives (for better 
readability, the points are just used by their number, e.g., Ff — F(a, vq + Avq, h)) 



F a = 



Ft, 



F n 



p — 



Fh„ — 



-3F\ + 4F 2 - F 3 

2Aa 
Fj - Fe 
2Av 

2Ah 

-3(-3Fi + 4F 2 - F 3 ) + 4(-3F 2 + 4F 3 - F 4 ) - (-3F 3 + 4F 4 - F 5 ) 

4(Act) 2 

(-3F 7 + 4^ - Fil) - (-3F 6 + 4^ - F 10 ) 



4AcrAw 
(-3Fi3 + 4F 15 - F 17 ) - (-3Fi2 
4AaAh 



4F 14 - F- 



16! 
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